Predicting polarization and nonlinear dielectric response of arbitrary perovskite 

superlattice sequences 
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We carry out first-principles calculations of the nonlinear dielectric response of short-period ferro- 
electric superlattices. We compute and store not only the total polarization, but also the Wannier- 
based polarizations of individual atomic layers, as a function of electric displacement field, and 
use this information to construct a model capable of predicting the nonlinear dielectric response of 
an arbitrary superlattice sequence. We demonstrate the successful application of our approach to 
superlattices composed of SrTiOa, CaTiOa, and BaTiOa layers. 

PACS numbers: 77.22.-d, 77.22.Ej, 77.80.-e, 77.84.Lf 



The development of advanced methods for layer-by- 
layer epitaxial growth of multicomponent perovskite su- 
perlattice structures has generated excitement both 
because of the intriguing materials physics that comes 
into play, and because of potential applications in non- 
volatile ferroelectric memories, piezoactuators and sen- 
sors, and magnetoelectric devices 0- To guide experi- 
mental exploration of this greatly expanded class of mate- 
rials, there is a critical need for atomic-scale understand- 
ing and modeling of key structural and functional prop- 
erties, particularly polarization and dielectric response. 

First-principles methods have allowed for the direct 
quantitative computation of such material-specific infor- 
mation for representative perovskite superlattices [^, |3| ■ 
However, such calculations are limited to relatively short- 
period superlattices, of the order of ten unit cells. First- 
principles modeling can extend our theoretical capabil- 
ity so that one can make predictions about arbitrary 
stacking sequences and elucidate the physics behind the 
novel behavior of superlattices. In particular, substantial 
progress has recently been made in isolating and study- 
ing the effects of the epitaxial strain on film structure, 
polarization, and piezoelectric properties (5|, Q . 

It is clear, however, that it is electrostatic effects that 
dominate the physics of superlattices built from ferro- 
electric (e.g., BaTiOs) and incipient ferroelectric (e.g., 
SrTiOs) constituents. In previous first-principles mod- 
els, these effects were included in an approximate way, 
either by describing the layers in terms of their bulk lin- 
ear dielectric properties (3|, or by imposing a constant- 
polarization layer-to-layer constraint that only roughly 
captures the effects of the internal electric fields . Fur- 
thermore, most previous first-principles calculations, us- 
ing the periodic boundary conditions implicit in ordinary 
implementations, give results only for zero applied elec- 
tric field. Since much of the interest in perovskite su- 
perlattices lies in their use in capacitor structures whose 
performance relies on their nonlinear dielectric behavior 
under bias voltage, a more fundamental methodology ca- 



pable of capturing such effects is urgently needed. 

In this Letter, we present a rigorous first-principles 
treatment allowing computation and modeling of the 
electrostatics and non-zero electric-field response of per- 
ovskite oxide superlattices. Our approach is based on 
a recently-developed Wannier-based formulation of the 
layer polarizations in perovskite superlattices 0] in com- 
bination with methods for treating insulators in finite 
electric fields 0, [13, El, [13. Crucially, we choose here to 



work at fixed electric displacement field [13J, and show 
that this gives a clean separation between long-range 
Coulomb interactions and short-range interfacial effects. 
As we demonstrate through application to superlattices 
composed of three ABO3 perovskite constituents, the re- 
sulting model yields, for arbitrary stacking sequences, 
quantitative predictions of polarization and nonlinear di- 
electric response with ab-initio accuracy, thus enabling 
the theory-driven search of the full range of superlattice 
sequences for novel or optimized properties. 

The construction of our model begins with the decom- 
position of the superlattice into atomic layers, specifi- 
cally into AO and BO2 layers alternating along [001]. 
The individual layer polarizations (LP) for each AO 
and BO2 layer j are computed using the Wannier-based 
method of Ref. [8], and recorded as functions Pj{D) of 
the displacement field D using a constrained-D first- 
principles implementation [l3|. This choice is appropri- 
ate because (i) D is constant throughout the superccU 
(V ■ D = 47r/9frco — 0), while the local macroscopic £ and 
P generally vary, and (ii) imposing constant-!? electrical 
boundary conditions has the virtue of making the force- 
constant matrix of the quasi-one-dimensional superlat- 
tice short-ranged in real space. This "locality principle" 
implies that one may expect the Pj{D) to depend only 
on the local compositional environment comprising the 
layer itself and few nearby neighbors. For any given su- 
perlattice, the total polarization, which is a quantity of 
central interest, is given by P{D) = c{D)^^ J2j Pji^)- I* 
is also straightforward to obtain the electric equation of 
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state in other forms, e.g., D{P) by numerical inversion, 
and £{P) = D{P) - AttP. 

We demonstrate the method through apphcation to su- 
perlattices comprised of an arbitrary sequence of SrTiOa, 
BaTiOa and CaTiOs layers, grown in the (001) di- 
rection with an in-plane lattice constant oq = 7.275 
bohr, our theoretical equilibrium lattice constant of bulk 
SrTiOs. We assume 1x1 in-plane periodicity and tetrag- 
onal Pimm symmetry, thus neglecting possible intermix- 
ing, nanodomain formation, or the appearance of antifer- 
roelectric or octahedral-tilting distortions. Some of these 
effects may be important in real superlattices and will 
deserve attention in future extensions, but our focus in 
the present work is to isolate and study the effects associ- 
ated with the complex stacking of ideal layers. Consistent 
with the PAmm symmetry, D is taken along the z axis, 
ranging in steps from —0.32 to 0.32 C/m^. 

First-principles calculations to optimize the struc- 
ture for fixed D [l3^ and to obtain the layer polar- 
izations Pj{D) and lattice constants c{D) were per- 
formed on a database of superlattices using the Lautrec 
code package, which implements plane-wave calcula- 
tions in the PAW framework 1^ in the LDA approx- 
imation [l5| . The polarization and its coupling to 
the electric field is handled by an efficient real-space 
Wannier formulation (12] . We used a plane- wave en- 
ergy cutoff of 40 Ry and a 6x6x2 Monkhorst-Pack k- 
mesh. The database of supcrlattice structures contains 
all one- and two-component period-4 supercells (BBBB, 
SSSS, CCCC, BBBS, BBSS, BSBS, BSSS, CCCS, CCSS, 
CSCS, CSSS, BBBC, BBCC, BCBC, BCCC), and one 
three-component supcrlattice SSBC, where C, S, and B 
refer to CaTiOs, SrTiOs, and BaTiOa layers respectively. 

Representative Pj{D) curves are presented in Fig. [1] 
It is striking that the LP curves separate, as expected 
from our locality principle, into clusters depending on 
the nearest-layer chemical environment (color-coded for 
comparison). However, the differences among the curves 
within a cluster are still too large to neglect, especially 
for Ti02 layers, indicating that an accurate model must 
include further-neighbor interactions as well. The effects 
of local inversion symmetry breaking are also clearly visi- 
ble. For example, a BaO layer in the middle of a CBB se- 
quence has a large and positive LP even at D=0. Smaller 
shifts arise from the second-neighbor environment, e.g., 
for the central Ti02 layer in a BBSS sequence. 

We now introduce a cluster expansion (CE) JJi] for the 
environment dependence of the p;({s}; D) as 

J'l , 5?) 
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where the J are D-dependent effective cluster interac- 



FIG. 1: (Color online.) Dependence of layer polarizations 
on chemical environment for (a) BaO planes (relative to a 
BaO plane in bulk BaTiOs), and (b) Ti02 planes (relative 
to an average of Ti02 planes in bulk BaTiOa and SrTiOs). 
C, S, B, and T refer to CaO, SrO, BaO, and Ti02 layers, 
respectively. First-principles results and model fittings are 
denoted by symbols and solid lines respectively. 



tions (ECI) to be determined from fitting to the first- 
principles database. We choose "spin" variables Si = —1, 
0, and 1 to identify AO layer i as CaO, SrO, and BaO 
repsectively, to reflect the fact that Sr is midway between 
Ca and Ba in the Periodic Table. Thus, insofar as Sr acts 
like an average of Ca and Ba atoms, each appearance of 
a squared spin variable in a term of the CE makes it 
more likely that the term can be neglected. 

We therefore approached the truncation and fitting of 
the model of Eq. ([T]) with three principles in mind: (1) 
the importance of an n-body term is expected to decrease 
rapidly with n; (2) the dependence of pi on Si should de- 
cay rapidly with the distance between layers I and i; and 
(3) coefficients with prime superscripts, corresponding to 
"higher- level" spin variables s^, should be less important 
than those without. 

Translation and spatial-inversion symmetry imply 

that Jl^l+rn{D) = ~Jl,l-m{ — D), Jl,l+mU+n{D) = 

— Ji,i-rn,i-n{—D), ctc, independent of I. It is there- 
fore natural to define Jm = {Ji,i—m + Ji,i+m)/'^ and 
J7m = {Ji,i~m — Ji,i+m) 1"^, and similarly for two-body and 
higher terms. Correspondingly, we separate Eq. H]) into 

parts p^p^Q (D) and PxToa (^) ^^^^ ^'^'^ ^^'^ parts 

Pao(^) ^^"^ PTi02(^) ^^'^ even in D, reflecting the 
inversion-symmetry-conserving and inversion-symmetry- 
breaking characters of the local environment respectively. 
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TABLE I: Fitted linear-in-D term of eflective cluster interac- 
tions for AO layers as defined in Eq. [2] 





ECI 


Value 


ECI 


Value 


Zero-body 


J 


2.2771 






One-body 


Jo 


0.1113 


Jo 


0.0819 




Ji 


0.0034 


Ji 


0.0007 






-0.0018 






Two-body 


Jq\ 


0.0197 


Jo2 


0.0031 




Jll 


0.0026 


J\1 


0.0013 



To give a sense of the form of the resulting expressions, 
the odd part for AO layers becomes 

pio^ = J + JoSo + Jo«o + J'l'lsi +sr) 

+ J2(S2 + S2) + Joi(siSo + SqSi) 
+ J02(S2S0 + S0S2) + JllSlSl 

+ Jll{s-2Sl^ SXS2) (2) 

where n=—n and layer is the one whose LP is being 
expanded. Similar expressions for p^O' P^.n ^ 1 and J3^02 
are given in the supplementary material [l7| . The super- 
cell lattice constant c{D) is correspondingly expanded as 

c = ^(Ci -I- CiSj + Czs] + C4SJSJ+1) , (3) 
3 

where Ci(£'), C2(-D), and Cj,(jy) assign to each layer its 
bulk c{D), and only the C^^D) term includes true su- 
perlattice effects. The J{D) and C{D) parameters are 
expressed as fifth and fourth-order Taylor expansions in 
D respectively, with the Taylor coefficients obtained by 
fitting to the first-principles calculations of Pj{D) and 
c{D) for superlattices in the database. 

The choice of terms to include in Eq. ([2]) and in the 

corresponding expressions for p^^Q , Pxi02 ' ^^'-^ PTi02 E3 
have been obtained using linear regression techniques to 
identify higher terms that could be deleted without sig- 
nificantly degrading the quality of the fit. The linear- 
in-Z? coefficients of the ECIs in Eq. ([2|) are presented in 
Table HI confirming our expectation that terms of higher 
body, longer range, and higher level tend to be less im- 
portant. Tables listing values of all of the ECI coefficients 
are provided in the supplementary material . 

The quality of the fit is excellent; the overall RMS 
error in pj {D) values relative to first-principles results is 
2 X 10~^'*C/m for structures in the database. This is 
illustrated in Fig. [1] where the solid lines representing 
the fitted functions can be seen to pass quite accurately 
through the first-principles symbols. The quality of the 
fit is similar for other cases, not plotted. 

The model can now be used to predict the nonlinear 
dielectric and piezoelectric properties of arbitrary super- 
lattice sequences. To illustrate the quality of the fit for 
supercell configurations that were not included in the fit, 
we compare in Fig.[2]the first-principles P{£) curves with 
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FIG. 2: (Color online.) Model prediction (solid lines) and 
first-principles calculations (symbols) of £ vs. P for 1S1B2C 
superlattice. 



the model fits for the tri-color 1S1B2C supercell. The 
P{£) curves are seen to be in excellent agreement. The 
arrows in Fig. [2] indicate £=0 solutions corresponding to 
P=—Ps in the preferred (stable) phase, -P=Punst at the 
unstable saddle point, and P=Pmeta in the metastable 
phase. Note that |Ps| ^ |Pmcta| and Punst be- 
cause the superlattice has broken inversion symmetry. 
The model predicts Pg, Punst, and Pmeta values of —0.19, 
0.04, and 0.16 C/m^, to be compared with direct first- 
principles values of —0.20, 0.04, and 0.16 C/m^, respec- 
tively. While the inversion-symmetry-breaking effects are 
subtle, they are critical for tuning certain ferroelectric 
properties [lil . [l9l | , and it is gratifying that they are ob- 
tained accurately by our model. 

To demonstrate the ability of our model to predict 
properties of long-period superlattices that would be 
impractical for direct first-principles calculations, we 
present our model predictions for the spontaneous polar- 
izations and dielectric responses (including the piezoelec- 
trically mediated component) of nSnBnC superlattices in 
Table |TT] and Fig. [3] respectively. Because of the broken 
inversion symmetry, one polarization direction is favored, 
with a different magnitude, over the other. The polariza- 
tions approach the bulk value of 0.274 C/m^ for large n, 
but with decreasing n we find a progressive suppression 
of the polarization and an enhancement of the asymme- 
try until, at n—1, the system becomes paraelectric, with 
a single minimum. These effects are also evident in the 
dielectric response curves shown in Fig. [31 where it is 



TABLE II: Predicted magnitude of polarization (C/m^) for 
preferred (Ps) and metastable (Pmcta) polarization states in 
nSnBnC superlattices. 





71=1 


n=2 


n=4 


n=8 


n=oo 


Ps 


-0.040 


-0.198 


-0.250 


-0.267 


-0.274 


-fmcta 




0.178 


0.237 


0.262 


0.274 
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FIG. 3: (Color online.) Dielectric response of nSnBnC super- 
lattice in (a) paraelectric and (b) ferroelectric regime. 

clear that the curves lack reflection symmetry about the 
vertical axis, and the system is seen to cross from ferro- 
electric to paraelectric behavior between n—2 and n=\. 

To gain more insight into interface effects in these su- 
perlattices, we can characterize each of the six kinds of 
interfaces by an interface dipole, extracted from our first- 
principles model as follows. Using S/B as an example, we 
imagine an infinite stack ...SSSBBB... with one interface, 
and define 

pUD) = Y.pm~pf\D) (4) 

i 

where i runs over AO and Ti02 layers, Pi{D) is the actual 
layer polarization predicted by our model, and p^^^ (D) is 
the polarization of that layer type in its own bulk environ- 
ment. (For the central Ti02 layer, p'^"^ (D) is the average 
of the bulk S and B values.) Because of the short-range 
nature of the model, the sum only needs to run over a 
few layers near the interface. 

The resulting Pint(-D) curves are presented in Fig. H) 
First, note that the curves tend to have a negative 
slope at D=0; this reflects the fact that the presence 
of interfaces tends to suppress the ferroelectricity, as 
was also evident in Table [Til and Fig. [31 Second, each 
pair such as BS and SB are related by the symmetry 
(-D) = —pf^{—D) required by the condition that the 
overall P{D) of a bicolor nSnB superlattice must be odd 
in D. Third, the inversion symmetry breaking, relevant 
to tricolor superlattices such as nBnSnC, is evident in 
the failure of the pf^tiD) + pf^{D) +p^^{D) curve to 
pass through the origin in the inset of Fig. [Jj 

The utility of the interface dipole concept is that, for 
any superlattice in which the interfaces are never sep- 
arated by less than three unit cells (as determined by 
the range of our model), P{D) can be calculated just 
by summing the bulk contribution for each layer and 
then adding the contribution from each interface. This 
prescription yields a simplified but rigorously equivalent 
model that can be used in such cases. Thus, for exam- 
ple, the inversion-symmetry-breaking effects in IBmSnC 
superlattices are captured exactly by the simplified pre- 
scription as long as m, n > 3. 
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FIG. 4: (Color online.) Model interface dipole densities. 

In summary, we have shown how a model can be ex- 
tracted from first-principles calculations on short-period 
superlattices and used to make quantitatively accurate 
predictions of nonlinear dielectric and piezoelectric re- 
sponses over a range of applied fields for arbitrary su- 
perlattice sequences. The treatment of electrostatic ef- 
fects is rigorous, a key aspect being the choice of the 
displacement field as the fundamental electrical variable 
so as to keep interlayer interactions short-ranged. The 
approach can be straightforwardly generalized to include 
dependence on epitaxial strain. Such an approach can 
play an important role in enabling the design of multi- 
functional ferroelectric superlattices with desired polar- 
ization, piezoelectric or dielectric responses. 
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